###############################################################################################
## Replication R Code for
## Babies across Borders: The Political Economy of International Child Adoption
## Authors: Asif Efrat, David Leblang, Steven Liao, Sonal S. Pandya
###############################################################################################

###############################################################################################
## Code date: January 24, 2015
## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Please contact Steven Liao <stevenliao@virginia.edu> for questions
################################################################################################

############
## Setup
############
# clear all objects in workspace
rm(list = ls())

# set scientific notation digits
options(scipen = 3)

# load packages
# sltool is a package in development by Steven Liao, 
# it includes various functions used in this study
# to install please follow instructions at <https://github.com/stevenliaotw/sltools>
pkg <- c("foreign", "plyr", "dplyr", "ggplot2", "maps", "geosphere", 
         "countrycode", "grid", "stringr", "rworldmap", "rworldxtra", 
         "RColorBrewer", "pscl", "sltools", "reporttools", "gridExtra",
         "texreg", "xtable")
lapply(pkg, require, character.only = TRUE)

# set directories
MAIN_DIR <- "/Users/stevenliao/Dropbox/international_adoption/final/replication"
# MAIN_DIR <- "replace above with your own directory here"
# create 2 folders named "data" and "out" within your main directory
# place replication data in the "data" folder
# the following code will load data from the "data" folder and output figures in the "out" folder
OUT_DIR <- paste(MAIN_DIR, "out", sep = "/")
DATA_DIR <- paste(MAIN_DIR, "data", sep = "/")

# load data
load(paste(DATA_DIR, "merge.RData", sep = "/"))
load(paste(DATA_DIR, "adopt.mi.RData", sep = "/"))
load(paste(DATA_DIR, "adopt.merge.RData", sep = "/"))

# set output directory
setwd(OUT_DIR)

# set seed
set.seed(1234)

# create MI data subsets excluding years 1991 and 1992
data.mi.m9192 <- llply(data.mi$imputations, function(x)
  return(filter(x, year >= 1993)))

# create MI data subsets excluding dyadic adoptions < 10
data.mi.m10 <- llply(data.mi.m9192, function(x)
  return(filter(x, adopt_tot >= 10)))

# subset raw data 1993-2010
data.9310 <- filter(data, year >= 1993 & year <= 2010)

###############################################################
## Fig 8. Global Count of International Adoptions, 1991-2010.
###############################################################
# get total adoptions by year
adopt.year <- adopt.merge %>% 
  group_by(year) %>%
  summarise(adopt_yearly = sum(adopt_tot, na.rm = TRUE)) %>%
  filter(year >= 1991 & year <= 2010)

# load library
pdf(file = "adoption.yearly.pdf", width = 6, height = 4)
adopt.yearly.f <- ggplot(adopt.year, aes(x = year, y = adopt_yearly)) +
  geom_line() +
  ggtitle("") +
  xlab("Year") + ylab("Number of Adoptions") +
  scale_x_continuous(breaks = c(1991, 1995, 2000, 2005, 2010)) +
  geom_vline(xintercept = 2004, colour = "red") + theme_bw()
adopt.yearly.f
dev.off()

########################################################################
## Fig 1. Number of Adoptees Sent to the United States, 1991-2010.
########################################################################
# some functions
checkDateLine <- function(l){
  n<-0
  k<-length(l)
  k<-k-1
  for (j in 1:k){
    n[j] <- l[j+1] - l[j]
  }
  n <- abs(n)
  m<-max(n, rm.na=TRUE)
  ifelse(m > 30, TRUE, FALSE)
}
clean.Inter <- function(p1, p2, n, addStartEnd){
  inter <- gcIntermediate(p1, p2, n = n, addStartEnd = addStartEnd)
  if (checkDateLine(inter[,1])){
    m1 <- midPoint(p1, p2)
    m1[,1] <- (m1[,1]+180)%%360 - 180
    a1 <- antipode(m1)
    l1 <- gcIntermediate(p1, a1, n = n, addStartEnd = addStartEnd)
    l2 <- gcIntermediate(a1, p2, n = n, addStartEnd = addStartEnd)
    l3 <- rbind(l1, l2)
    l3
  }
  else{
    inter
  }
}

# subset
adopt.usa <- filter(adopt.merge, iso3c_d == "USA")

# aggregate dyad adopt_tot 1991-2010
adopt.usa.s <- adopt.usa %>%
  group_by(iso3c_o, iso3c_d) %>%
  summarise(adopt_dyad_us = sum(adopt_tot, na.rm = TRUE))

# load coordinates data
coord <- read.dta(paste(DATA_DIR, "dyad.map.dta", sep = "/"))

# subset
coord <- coord %>%
  dplyr::select(iso3_d, iso3_o, longitude_o, latitude_o, longitude_d, latitude_d) %>%
  group_by(iso3_d, iso3_o) %>%
  slice(1)

# merge with US adoption data
adopt.usa.s <- merge(adopt.usa.s, coord, by.x = c("iso3c_d", "iso3c_o"), by.y = c("iso3_d", "iso3_o"), all.x = TRUE)

# exclude entities with no coords
filter(adopt.usa.s, is.na(longitude_o))
adopt.usa.s <- filter(adopt.usa.s, !is.na(longitude_o))

# set colors and line width
#pal <- colorRampPalette(c("#f2f2f2", "black"))
#pal <- colorRampPalette(c("#f2f2f2", "red"))
#pal <- colorRampPalette(c("#f2f2f2", "blue"))
#pal <- colorRampPalette(c("#f2f2f2"))
pal <- colorRampPalette(c("red"))
#pal <- colorRampPalette(c("#00FF00", "#FF0000"))
colors <- pal(100)
inter.lwd <- seq(0, 4, length = 100)

# plot map
pdf(file = "flowmap_USA.pdf", width = 7, height = 4.2)
par(mai = c(0, 0, 0, 0))
map("world", col = "#f2f2f2", fill = TRUE, bg = "white", lwd = 0.2)
#map("world", col="#191919", fill=TRUE, bg="#736F6E", lwd=0.05, mai = c(0, 0, 0.4, 0))
#map("world", col="black", fill=FALSE, bg="white", lwd=0.15, mar = c(0, 0, 0.2, 0))

# sort
adopt.usa.s <- arrange(adopt.usa.s, adopt_dyad_us)

# extract max count for scaling line width and color
maxcnt <- max(adopt.usa.s$adopt_dyad_us, na.rm = TRUE)

# draw lines for each dyad value
for (j in 1:length(adopt.usa.s$iso3c_d)) {
  p1 <- c(adopt.usa.s[j,]$longitude_o, adopt.usa.s[j,]$latitude_o)
  p2 <- c(adopt.usa.s[j,]$longitude_d, adopt.usa.s[j,]$latitude_d)
  inter <- clean.Inter(p1, p2, n = 800, addStartEnd=TRUE)
  colindex <- round((adopt.usa.s[j,]$adopt_dyad_us/maxcnt)*length(colors))
  lwdindex <- round((adopt.usa.s[j,]$adopt_dyad_us/maxcnt)*length(inter.lwd))
  lines(inter, col=colors[colindex], lwd = inter.lwd[lwdindex])
}
box()
mtext("", col = "Black", cex = 0.7, side = 3, line = 0.7)
segments(-180, -40, -170, -40, 
         lwd = inter.lwd[round((filter(adopt.usa.s, iso3c_o == "CHN")$adopt_dyad_us/maxcnt)*length(colors))],
         col = colors[round((filter(adopt.usa.s, iso3c_o == "CHN")$adopt_dyad_us/maxcnt)*length(inter.lwd))])
text(-170, -40, paste(filter(adopt.usa.s, iso3c_o == "CHN")$adopt_dyad_us, "(China)"), pos = 4, cex = 0.6)
segments(-180, -50, -170, -50, 
         lwd = inter.lwd[round((filter(adopt.usa.s, iso3c_o == "GTM")$adopt_dyad_us/maxcnt)*length(colors))],
         col = colors[round((filter(adopt.usa.s, iso3c_o == "GTM")$adopt_dyad_us/maxcnt)*length(inter.lwd))])
text(-170, -50, paste(filter(adopt.usa.s, iso3c_o == "GTM")$adopt_dyad_us, "(Guatemala)"), pos = 4, cex = 0.6)
segments(-180, -60, -170, -60, 
         lwd = inter.lwd[round((filter(adopt.usa.s, iso3c_o == "MEX")$adopt_dyad_us/maxcnt)*length(colors))],
         col = colors[round((filter(adopt.usa.s, iso3c_o == "MEX")$adopt_dyad_us/maxcnt)*length(inter.lwd))])
text(-170, -60, paste(filter(adopt.usa.s, iso3c_o == "MEX")$adopt_dyad_us, "(Mexico)"), pos = 4,cex = 0.6)
dev.off()

########################################################################
##Fig 2. Number of Adoptees by Receiving Country, 1991-2010.
########################################################################
# extract yearly adoption data for destination country, 1991-2010
adopt.d.panel <- adopt.merge %>%
  group_by(iso3c_d, year) %>%
  summarise(adopt_d_yearly = sum(adopt_tot, na.rm = TRUE)) %>%
  filter(year >= 1991 & year <= 2010)

# add country names
adopt.d.panel$cty_d <- countrycode(adopt.d.panel$iso3c_d, "iso3c", "country.name")

# plot
adopt.d.yearly.f <- ggplot(adopt.d.panel, aes(x = year, y = adopt_d_yearly)) +
  geom_line() +
  ggtitle("") +
  xlab("Year") + ylab("Adoptees Received") +
  scale_x_continuous(breaks = c(1991, 1995, 2000, 2005, 2010)) +
  facet_wrap(~cty_d, ncol = 5) + 
  theme_bw() +
  theme(panel.margin = unit(0.5, "lines"))
pdf(file = "adoption.d.yearly.pdf", width = 13, height = 9)
adopt.d.yearly.f
dev.off()

########################################################################
## Fig 3. Number of Adoptees from 35 Select Sending Countries, 1991-2010.
########################################################################
# extract yearly adoption data for origin country, 1991-2010
adopt.o.panel <- adopt.merge %>%
  group_by(iso3c_o, year) %>%
  summarise(adopt_o_yearly = sum(adopt_tot, na.rm = TRUE)) %>%
  filter(year >= 1991 & year <= 2010)

# add country name
adopt.o.panel$cty_o <- countrycode(adopt.o.panel$iso3c_o, "iso3c", "country.name")

# plot panel data for all sending countries
adopt.o.yearly.f <- ggplot(adopt.o.panel, aes(x = year, y = adopt_o_yearly)) +
  geom_line() +
  ggtitle("") +
  xlab("Year") + ylab("Adoptees Sent") +
  scale_x_continuous(breaks = c(1991, 1995, 2000, 2005, 2010)) +
  facet_wrap(~cty_o, ncol = 5) + 
  theme_bw() + 
  theme(panel.margin = unit(0.5, "lines"))
pdf(file = "adoption.o.yearly.pdf", width = 20, height = 100)
adopt.o.yearly.f
dev.off()

# calculate total adoptees sent by origin
adopt.o.panel <- mutate(adopt.o.panel, adopt_o_tot = sum(adopt_o_yearly, na.rm = TRUE))

# correct names
adopt.o.panel$cty_o <- str_replace(adopt.o.panel$cty_o, "Taiwan, Province Of China", "Taiwan")
adopt.o.panel$cty_o <- str_replace(adopt.o.panel$cty_o, "Bolivia, Plurinational State Of", "Bolivia")
adopt.o.panel$cty_o <- str_replace(adopt.o.panel$cty_o, "Korea, Republic Of", "South Korea")

# subset countries sending 1500 or more adoptees
adopt.o.panel.sub <- filter(adopt.o.panel, adopt_o_tot >= 1500)

# plot
adopt.o.yearly.f <- ggplot(adopt.o.panel.sub, aes(x = year, y = adopt_o_yearly)) +
  geom_line() +
  ggtitle("") +
  xlab("Year") + ylab("Adoptees Sent") +
  scale_x_continuous(breaks = c(1991, 1995, 2000, 2005, 2010)) +
  facet_wrap(~cty_o, ncol = 5) + 
  theme_bw() + 
  theme(panel.margin = unit(0.5, "lines"))
pdf(file = "adoption.o.yearly.1500.pdf", width = 12, height = 15)
adopt.o.yearly.f
dev.off()

########################################################################
## Fig 4. Total Outflows and Inflows of Adoptees, 1991-2010
########################################################################
# generating monadic data and subset year
adopt.o.map <- adopt.merge %>%
  group_by(iso3c_o, year) %>%
  summarise(adopt_o = sum(adopt_tot, na.rm = TRUE),
            cty_o = head(cty_o, 1)) %>%
  filter(year >= 1991 & year <= 2010)

adopt.d.map <- adopt.merge %>%
  group_by(iso3c_d, year) %>%
  summarise(adopt_d = sum(adopt_tot, na.rm = TRUE),
            cty_d = head(cty_d, 1)) %>%
  filter(year >= 1991 & year <= 2010)

# aggregate all adoptions per origin countries, 1991-2010
adopt.o.map.1991.2010.s <- summarise(adopt.o.map,
                                     adopt_o_91_10 = sum(adopt_o, na.rm = TRUE),
                                     cty_o = head(cty_o, 1))

# aggregate all adoptions per destination countries, 1991-2010
adopt.d.map.1991.2010.s <- summarise(adopt.d.map,
                                     adopt_d_91_10 = sum(adopt_d, na.rm = TRUE),
                                     cty_d = head(cty_d, 1))

# merge adoption data with map data
# origin countries
adopt.o.91.10 <- joinCountryData2Map(adopt.o.map.1991.2010.s, 
                                     joinCode = "ISO3", 
                                     nameJoinColumn = "iso3c_o",
                                     verbose = TRUE)

# destination countries
adopt.d.91.10 <- joinCountryData2Map(adopt.d.map.1991.2010.s, 
                                     joinCode = "ISO3", 
                                     nameJoinColumn = "iso3c_d",
                                     verbose = TRUE)

# set bins for choropleth maps
adopt.breaks.snapshot.o <- c(0, 1, 10, 100, 1000, 10000, 100000, 150000)
adopt.breaks.snapshot.d <- c(0, 1, 10, 100, 1000, 10000, 100000, 500000)

# plot sending countries
mapDevice(device = "pdf", width = 7, 
          file = "choropleth_snapshot.pdf",
          rows = 2, columns = 1, 
          mai = c(0.1, 0.1, 0.5, 1))
adopt.o.rworldmap.91.10 <- mapCountryData(adopt.o.91.10,
                                          nameColumnToPlot = "adopt_o_91_10",
                                          mapRegion = "world",
                                          numCats = 7,
                                          catMethod = adopt.breaks.snapshot.o,
                                          missingCountryCol = "white",
                                          mapTitle = "Sending Countries' Aggregate Adoption Outflows, 1991-2010",
                                          addLegend = FALSE,
                                          colourPalette = brewer.pal(7, "Blues"))
do.call(addMapLegend, c(adopt.o.rworldmap.91.10,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 5,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
adopt.d.rworldmap.91.10 <- mapCountryData(adopt.d.91.10,
                                          nameColumnToPlot = "adopt_d_91_10",
                                          mapRegion = "world",
                                          numCats = 7,
                                          catMethod = adopt.breaks.snapshot.d,
                                          missingCountryCol = "white",
                                          mapTitle = "Receiving Countries' Aggregate Adoption Inflows, 1991-2010",
                                          addLegend = FALSE,
                                          colourPalette = brewer.pal(7, "Blues"))
do.call(addMapLegend, c(adopt.d.rworldmap.91.10,
                        legendLabels = "all",
                        labelFontSize = 0.8,
                        legendMar = 5,
                        legendWidth = 0.5,
                        legendShrink = 0.9,
                        horizontal = FALSE))
dev.off()

#####################
## Fit Full Models
#####################
# make vector of covariates
covariate.full <- c(## transaction costs unique to adoption
                    "reg_quality_o_L1",
                    "execnat_o_L1",
                    "log(adopt_o_cumsum_L1 + 1)",
                    "hague_adopt_entry_both_L1",

                    ## transaction costs general to migration
                    "comlang_ethno",
                    "log(Istock_L1 + 1)",
                    "colony",
                    "log(distcap)",
                    
                    ## supply of children (controls)
                    "log(pop14_mil_L1)",
                    "log(rgdppc_exp_L1)",
                    "immune_measles_L1",
                    "conflict_o_L1",
                    "islam_law_o",
                    
                    ## destination country fixed effects
                    "factor(ifs_d)", "factor(year)")

# make formula for full model
formula.full <- formula(paste("adopt_tot ~ ", paste(covariate.full, collapse = " + ")))
formula.full

## fit full hurdle model (logit + negative binomial) to 10 MI datasets, period 1993-2010
hurdle.full.m9192.out <- llply(data.mi.m9192, function(x)
  return(hurdle(formula.full,
                dist = "negbin",
                zero.dist = "binomial", 
                link = "logit",
                data = x)))

# Get MI results
hurdle.full.m9192.results <- combineHurdleMI(fitted.obj = hurdle.full.m9192.out)
hurdle.full.m9192.results

########################################################################
## Fig 5. Hurdle Model of the Determinants of International Adoption.
########################################################################
# some functions to prepare parameters for plotting
prepFigZeroPars <- function(hurdle.out, alpha){
  # extract coefficients
  coef <- hurdle.out$zero.coef
  # set multiplier based on alpha and degrees of freedom
  multiplier <- qt(1 - alpha/2, df = hurdle.out$dof.zero)
  # calculate high and low values
  coef.high <- coef + hurdle.out$zero.se*multiplier
  coef.low <- coef - hurdle.out$zero.se*multiplier
  # omit factors
  coef <- coef[!grepl("factor", names(coef))]
  coef.high <- coef.high[!grepl("factor", names(coef.high))]
  coef.low <- coef.low[!grepl("factor", names(coef.low))]
  # omit intercept
  coef <- coef[!grepl("Intercept", names(coef))]
  coef.high <- coef.high[!grepl("Intercept", names(coef.high))]
  coef.low <- coef.low[!grepl("Intercept", names(coef.low))]
  # make list
  list("mean" = coef,
       "high" = coef.high,
       "low" = coef.low)
}

prepFigCountPars <- function(hurdle.out, alpha){
  # extract coefficients
  coef <- hurdle.out$count.coef
  # set multiplier based on alpha and degrees of freedom
  multiplier <- qt(1 - alpha/2, df = hurdle.out$dof.count)
  # calculate high and low values
  coef.high <- coef + hurdle.out$count.se*multiplier
  coef.low <- coef - hurdle.out$count.se*multiplier
  # omit factors
  coef <- coef[!grepl("factor", names(coef))]
  coef.high <- coef.high[!grepl("factor", names(coef.high))]
  coef.low <- coef.low[!grepl("factor", names(coef.low))]
  # omit intercept
  coef <- coef[!grepl("Intercept", names(coef))]
  coef.high <- coef.high[!grepl("Intercept", names(coef.high))]
  coef.low <- coef.low[!grepl("Intercept", names(coef.low))]
  # make list
  list("mean" = coef,
       "high" = coef.high,
       "low" = coef.low)
}

# extract cleaned parameters
coef.zero.m9192.fig <- prepFigZeroPars(hurdle.out = hurdle.full.m9192.results, 
                                       alpha = 0.05)
coef.count.m9192.fig <- prepFigCountPars(hurdle.out = hurdle.full.m9192.results, 
                                         alpha = 0.05)

# assign full variable names
coef.names.m9192 <- c("Requlatory Quality",
                      "Nationalism",
                      "ln(Cumulative Adoption)",
                      "Hague Adoption (both)",
                      "Language Commonality",
                      "ln(Migrant Stock)",
                      "Colonial Ties",
                      "ln(Distance)",
                      "ln(Youth Population)",
                      "ln(Real GDP per Capita)",
                      "Measles Immunization",
                      "Armed Conflict",
                      "Islamic Law")

# prepare parameters for plot
# create indicator for y.axes, descending so that R orders vars from top to bottom on y-axis
y.axis.zero.m9192.coef <- c(length(coef.zero.m9192.fig$mean):1) 
y.axis.count.m9192.coef <- c(length(coef.count.m9192.fig$mean):1) 

# open pdf device
pdf("coefPlotCombine9310.pdf", height = 5, width = 9) 
# set plot panels
nPanels <- layout(rbind(c(1, 2)), heights = c(3), respect = FALSE)
# set margins for plot, leaving lots of room on left-margin (2nd number in margin command) for variable names
par(mar = c(2, 9, 2, 0)) 
# plot logit component coefs of hurdle model
plot(coef.zero.m9192.fig$mean, 
     y.axis.zero.m9192.coef, 
     # plot coefficients as points, turning off axes and labels.
     type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black", 
     # set limits of x-axis so that they include mins and maxs of coefficients + .95% confidence intervals
     xlim = c(round(min(coef.zero.m9192.fig$low), digits = 2), 
              round(max(coef.zero.m9192.fig$high), digits = 2)),
     ylim = c(min(y.axis.zero.m9192.coef), 
              max(y.axis.zero.m9192.coef)),
     xaxs = "r", main = "") 
# add segment for coef intervals
segments(coef.zero.m9192.fig$low, 
         y.axis.zero.m9192.coef, 
         coef.zero.m9192.fig$high, 
         y.axis.zero.m9192.coef, 
         lwd =  1.5, col="black")
# add axes
axis(1, at = seq(round(min(coef.zero.m9192.fig$low), digits = 2) - 0.03,
                 round(max(coef.zero.m9192.fig$high), digits = 2) + 0.13,
                 by = 0.1),
     labels =  seq(round(min(coef.zero.m9192.fig$low), digits = 2) - 0.03,
                   round(max(coef.zero.m9192.fig$high), digits = 2) + 0.13,
                   by = 0.1),
     # draw x-axis and labels with tick marks
     tick = TRUE, 
     cex.axis = .8, 
     # reduce label size, moves labels closer to tick marks
     mgp = c(2, .5, 0))
axis(2, at = y.axis.zero.m9192.coef, 
     label = coef.names.m9192, 
     las = 1, 
     tick = TRUE,
     cex.axis = .8)
# draw dotted line through 0
segments(0, 0, 0, 17, lty = 2) 
# annotate
mtext("Prob(Adoption)", side = 3)

# plot count component coefs of hurdle model
par(mar = c(2, 1, 2, 0.5))
plot(coef.count.m9192.fig$mean, 
     y.axis.count.m9192.coef, 
     type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black",
     xlim = c(floor(min(coef.count.m9192.fig$low)), 
              ceiling(max(coef.count.m9192.fig$high))),
     ylim = c(min(y.axis.count.m9192.coef), 
              max(y.axis.count.m9192.coef)),
     xaxs = "r", main = "") 
segments(coef.count.m9192.fig$low, 
         y.axis.count.m9192.coef, 
         coef.count.m9192.fig$high, 
         y.axis.count.m9192.coef, 
         lwd =  1.5, col="black") 
axis(1, at = seq(floor(min(coef.count.m9192.fig$low)), 
                 ceiling(max(coef.count.m9192.fig$high)), 
                 by = 0.1),
     labels =  seq(floor(min(coef.count.m9192.fig$low)), 
                   ceiling(max(coef.count.m9192.fig$high)), 
                   by = 0.1),
     tick = TRUE, 
     cex.axis = .8, 
     mgp = c(2,.5,0)) 
segments(0, 0 , 0, 17, lty = 2)
mtext("Adoption Count", side = 3)
text(0.51, 3, paste("n = ", hurdle.full.m9192.results$n), cex = 0.8)
text(0.65, 2, paste("logLik = ", round(hurdle.full.m9192.results[[9]], digits = 2)), cex = 0.8)
text(0.6, 1, paste("AIC = ", round(hurdle.full.m9192.results[[10]], digits = 2)), cex = 0.8)
dev.off()

################################################################################################################################################
## simulating effects based on Zelig method and hurdle model results for period 1993-2010
################################################################################################################################################
# define the values at which the variables are held constant, and the variable of interest that is allowed to vary
# bothhague_L1 = 0 to 1, and mean values for all other continuous variables, mode for dichotomous variables
# set default values
set.bothhague <- data.frame(intercept = 1,
                            ## transaction costs unique to adoption
                            reg_quality_o_L1 = mean(ldply(data.mi.m9192, function(x) mean(x$reg_quality_o_L1))[,2]),
                            execnat_o_L1 = 0,
                            "log(adopt_o_cumsum_L1 + 1)" = mean(ldply(data.mi.m9192, function(x) mean(log(x$adopt_o_cumsum_L1 + 1)))[,2]),
                            hague_adopt_entry_both_L1 = 0:1,
                            
                            ## transaction costs general to migration
                            comlang_ethno = 0,
                            "log(Istock_L1 + 1)" = mean(ldply(data.mi.m9192, function(x) mean(log(x$Istock_L1 + 1)))[,2]),
                            colony = 0,
                            "log(distcap)" = mean(ldply(data.mi.m9192, function(x) mean(log(x$distcap)))[,2]),
                            
                            ## controls
                            "log(pop14_mil_L1)" = mean(ldply(data.mi.m9192, function(x) mean(log(x$pop14_mil_L1)))[,2]),
                            "log(rgdppc_exp_L1)" = mean(ldply(data.mi.m9192, function(x) mean(log(x$rgdppc_exp_L1)))[,2]),
                            immune_measles_L1 = mean(ldply(data.mi.m9192, function(x) mean(x$immune_measles_L1))[,2]),
                            conflict_o_L1 = 0,
                            islam_law_o = 0)
set.bothhague[c("factor(ifs_d)112", "factor(ifs_d)124", "factor(ifs_d)128", "factor(ifs_d)132", "factor(ifs_d)134", "factor(ifs_d)136",
                "factor(ifs_d)138", "factor(ifs_d)142", "factor(ifs_d)144", "factor(ifs_d)146", "factor(ifs_d)156", "factor(ifs_d)172",
                "factor(ifs_d)176", "factor(ifs_d)178", "factor(ifs_d)184", "factor(ifs_d)193", "factor(ifs_d)196", "factor(ifs_d)436",
                #"factor(year)1992", "factor(year)1993", 
                "factor(year)1994", "factor(year)1995", "factor(year)1996", "factor(year)1997",
                "factor(year)1998", "factor(year)1999", "factor(year)2000", "factor(year)2001", "factor(year)2002", "factor(year)2003",
                "factor(year)2004", "factor(year)2005", "factor(year)2006", "factor(year)2007", "factor(year)2008", "factor(year)2009",
                "factor(year)2010")] <- 0
set.bothhague["factor(year)2000"] <- 1
dim(set.bothhague)

# simulate both hague effects
s.out.bothhague <- simHurdle(data = data.mi.m9192, 
                             x = set.bothhague[1,], 
                             x1 = set.bothhague[2,], 
                             hurdle.out = hurdle.full.m9192.out, 
                             num = 10000)
s.out.bothhague$fd.zero
s.out.bothhague$fd.count

# simulate reg_quality_o_L1
# set values
set.rq <- set.bothhague
# reset both hague variable to 0s
set.rq$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
mu.rq <- mean(ldply(data.mi.m9192, function(x) mean(x$reg_quality_o_L1))[,2])
sd.rq <- mean(ldply(data.mi.m9192, function(x) sd(x$reg_quality_o_L1))[,2])
set.rq$reg_quality_o_L1 <- c(mu.rq - sd.rq, mu.rq + sd.rq)
set.rq
# simulate
s.out.rq <- simHurdle(data = data.mi.m9192, 
                      x = set.rq[1,], 
                      x1 = set.rq[2,], 
                      hurdle.out = hurdle.full.m9192.out, 
                      num = 10000)
s.out.rq$fd.zero
s.out.rq$fd.count

# simulate Nationalism
# set values
set.nationalism <- set.bothhague
# reset both hague variable
set.nationalism$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
set.nationalism$execnat_o_L1 <- c(0, 1)
set.nationalism
# simulate
s.out.nationalism <- simHurdle(data = data.mi.m9192, 
                               x = set.nationalism[1,], 
                               x1 = set.nationalism[2,], 
                               hurdle.out = hurdle.full.m9192.out, 
                               num = 10000)
s.out.nationalism$fd.zero
s.out.nationalism$fd.count

# simulate Cumulative Adoption
# set values
set.ac <- set.bothhague
# reset both hague variable
set.ac$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
mu.ac <- mean(ldply(data.mi.m9192, function(x) mean(log(x$adopt_o_cumsum_L1 + 1)))[,2])
sd.ac <- mean(ldply(data.mi.m9192, function(x) sd(log(x$adopt_o_cumsum_L1 + 1)))[,2])
set.ac[, grep("*cumsum", names(set.ac))] <- c(mu.ac - sd.ac, mu.ac + sd.ac)
set.ac
# simulate
s.out.ac <- simHurdle(data = data.mi.m9192, 
                      x = set.ac[1,], 
                      x1 = set.ac[2,], 
                      hurdle.out = hurdle.full.m9192.out, 
                      num = 10000)
s.out.ac$fd.zero
s.out.ac$fd.count

# simulate Language Commonality (comlang_ethno_L1)
# set values
set.lc <- set.bothhague
# reset both hague variable
set.lc$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
set.lc$comlang_ethno <- c(0, 1)
set.lc
# simulate
s.out.lc <- simHurdle(data = data.mi.m9192, 
                      x = set.lc[1,], 
                      x1 = set.lc[2,], 
                      hurdle.out = hurdle.full.m9192.out, 
                      num = 10000)
s.out.lc$fd.zero
s.out.lc$fd.count

# simulate Migrant Stock (Istock_L1)
# set values
set.ms <- set.bothhague
# reset both hague variable
set.ms$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
mu.ms <- mean(ldply(data.mi.m9192, function(x) mean(log(x$Istock_L1 + 1)))[,2])
sd.ms <- mean(ldply(data.mi.m9192, function(x) sd(log(x$Istock_L1 + 1)))[,2])
set.ms[, grep("Istock", names(set.ms))] <- c(mu.ms - sd.ms, mu.ms + sd.ms)
set.ms
# simulate
s.out.ms <- simHurdle(data = data.mi.m9192, 
                      x = set.ms[1,], 
                      x1 = set.ms[2,], 
                      hurdle.out = hurdle.full.m9192.out, 
                      num = 10000)
s.out.ms$fd.zero
s.out.ms$fd.count

# simulate Colonial Ties (colony)
# set values
set.ct <- set.bothhague
# reset both hague variable
set.ct$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
set.ct$colony <- c(0, 1)
set.ct
# simulate
s.out.ct <- simHurdle(data = data.mi.m9192, 
                      x = set.ct[1,], 
                      x1 = set.ct[2,], 
                      hurdle.out = hurdle.full.m9192.out, 
                      num = 10000)
s.out.ct$fd.zero
s.out.ct$fd.count

# simulate Distance (distcap)
# set values
set.dist <- set.bothhague
# reset both hague variable
set.dist$hague_adopt_entry_both_L1 <- c(0, 0)
# set new variable of interest
mu.dist <- mean(ldply(data.mi.m9192, function(x) mean(log(x$distcap + 1)))[,2])
sd.dist <- mean(ldply(data.mi.m9192, function(x) sd(log(x$distcap + 1)))[,2])
set.dist[, grep("distcap", names(set.dist))] <- c(mu.dist - sd.dist, mu.dist + sd.dist)
set.dist
# simulate
s.out.dist <- simHurdle(data = data.mi.m9192, 
                        x = set.dist[1,], 
                        x1 = set.dist[2,], 
                        hurdle.out = hurdle.full.m9192.out, 
                        num = 10000)
s.out.dist$fd.zero
s.out.dist$fd.count

# combine all simulated first difference results
effect.negbin.data <- rbind(s.out.rq$fd.count[, 3],
                            s.out.nationalism$fd.count[, 3],
                            s.out.ac$fd.count[, 3],
                            s.out.bothhague$fd.count[, 3],
                            s.out.lc$fd.count[, 3],
                            s.out.ms$fd.count[, 3],
                            s.out.ct$fd.count[, 3],
                            s.out.dist$fd.count[, 3])
effect.zero.data <- rbind(s.out.rq$fd.zero[, 3],
                          s.out.nationalism$fd.zero[, 3],
                          s.out.ac$fd.zero[, 3],
                          s.out.bothhague$fd.zero[, 3],
                          s.out.lc$fd.zero[, 3],
                          s.out.ms$fd.zero[, 3],
                          s.out.ct$fd.zero[, 3],
                          s.out.dist$fd.zero[, 3])

# extract vectors of interest for effects plot
effect.negbin.vec <- effect.negbin.data[, 2]
effect.negbin.low.vec <- effect.negbin.data[, 1]
effect.negbin.high.vec <- effect.negbin.data[, 3]
effect.zero.vec <- effect.zero.data[, 2]
effect.zero.low.vec <- effect.zero.data[, 1]
effect.zero.high.vec <- effect.zero.data[, 3]

# assign full var names
effect.names <- c("Requlatory Quality",
                  "Nationalism",
                  "ln(Cumulative Adoption)",
                  "Both Hague",
                  "Language Commonality",
                  "ln(Migrant Stock)",
                  "Colonial Ties",
                  "ln(Distance)")

# prepare parameters for plot
y.axis.negbin.effect <- c(length(effect.negbin.vec):1)
y.axis.zero.effect <- c(length(effect.zero.vec):1) 

pdf("effectPlotCombine9310.pdf", height = 5, width = 9) 
nPanels <- layout(rbind(c(1, 2)), heights = c(3), respect = FALSE)
par(mar = c(2, 9, 2, 0)) 
plot(effect.zero.vec, 
     y.axis.zero.effect, 
     type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black", 
     xlim = c(round(min(effect.zero.low.vec), digits = 2), 
              round(max(effect.zero.high.vec), digits = 2)),
     ylim = c(min(y.axis.zero.effect), 
              max(y.axis.zero.effect)),
     xaxs = "r", main = "") 
segments(effect.zero.low.vec, 
         y.axis.zero.effect, 
         effect.zero.high.vec, 
         y.axis.zero.effect, 
         lwd =  1.5, 
         col="black") 
axis(1, at = seq(round(min(effect.zero.low.vec), 
                       digits = 2) - 0.04,
                 round(max(effect.zero.high.vec), 
                       digits = 2) + 0.01,
                 by = 0.15),
     labels =  seq(round(min(effect.zero.low.vec), 
                         digits = 2) - 0.04,
                   round(max(effect.zero.high.vec), 
                         digits = 2) + 0.01,
                   by = 0.15),
     tick = TRUE,
     cex.axis = .8, 
     mgp = c(2, .5, 0)) 
axis(2, at = y.axis.zero.effect, 
     label = effect.names, 
     las = 1, 
     tick = TRUE,
     cex.axis = .8) 
segments(0, 0, 0, 17, lty = 2)
mtext("Prob(Adoption)", side = 3)

par(mar = c(2, 1, 2, 0.5)) 
plot(effect.negbin.vec, 
     y.axis.negbin.effect, 
     type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black",
     xlim = c(floor(min(effect.negbin.low.vec)), 
              ceiling(max(effect.negbin.high.vec)) + 3),
     ylim = c(min(y.axis.negbin.effect), 
              max(y.axis.negbin.effect)),
     xaxs = "r", main = "") 
segments(effect.negbin.low.vec, 
         y.axis.negbin.effect, 
         effect.negbin.high.vec, 
         y.axis.negbin.effect, 
         lwd =  1.5, 
         col = "black")
axis(1, at = seq(floor(min(effect.negbin.low.vec)) - 1, 
                 ceiling(max(effect.negbin.high.vec)) + 2, by = 5),
     labels =  seq(floor(min(effect.negbin.low.vec)) - 1, 
                   ceiling(max(effect.negbin.high.vec)) + 2, by = 5),
     tick = TRUE, 
     cex.axis = .8, 
     mgp = c(2, .5, 0)) 
segments(0, 0 , 0, 17, lty = 2)
mtext("Adoption Count", side = 3)
dev.off()

##########################################
## Robustness Checks: Alternative Models
##########################################
# poisson
pois.full <- llply(data.mi.m9192, function(x)
  return(glm(formula.full,
             family = poisson,
             data = x)))
pois.full.results <- combineGlmMI(pois.full)
pois.full.results

# negative binomial
nb.full <- llply(data.mi.m9192, function(x)
  return(glm.nb(formula.full,
                data = x,
                control = glm.control(epsilon = 1e-8, maxit = 200, trace = FALSE))))
nb.full.results <- combineGlmMI(nb.full)
nb.full.results

# zero-inflated model
zi.full <- llply(data.mi.m9192, function(x)
  return(zeroinfl(formula.full,
                  dist = "negbin",
                  link = "logit",
                  data = x)))
zi.full.results <- combineZIcountMI(fitted.obj = zi.full)
zi.full.results

##########################################
## Robustness Checks: Alternative Indicators
##########################################
# control of corruption specification
covariate.corrupt <- c(## transaction costs unique to adoption
                        "control_corrupt_o_L1",
                        "execnat_o_L1",
                        "log(adopt_o_cumsum_L1 + 1)",
                        "hague_adopt_entry_both_L1",
                        
                        ## transaction costs general to migration
                        "comlang_ethno",
                        "log(Istock_L1 + 1)",
                        "colony",
                        "log(distcap)",
                        
                        ## supply of children (controls)
                        "log(pop14_mil_L1)",
                        "log(rgdppc_exp_L1)",
                        "immune_measles_L1",
                        "conflict_o_L1",
                        "islam_law_o",
                        
                        ## destination country fixed effects
                        "factor(ifs_d)", "factor(year)")
formula.corrupt <- formula(paste("adopt_tot ~ ", paste(covariate.corrupt, collapse = " + ")))
formula.corrupt

# jus soli specification
covariate.soli <- c(## transaction costs unique to adoption
  "reg_quality_o_L1",
  "Mjussolis_o_L1",
  "log(adopt_o_cumsum_L1 + 1)",
  "hague_adopt_entry_both_L1",

  ## transaction costs general to migration
  "comlang_ethno",
  "log(Istock_L1 + 1)",
  "colony",
  "log(distcap)",
  
  ## supply of children (controls)
  "log(pop14_mil_L1)",
  "log(rgdppc_exp_L1)",
  "immune_measles_L1",
  "conflict_o_L1",
  "islam_law_o",

  ## destination country fixed effects
  "factor(ifs_d)", "factor(year)")
formula.soli <- formula(paste("adopt_tot ~ ", paste(covariate.soli, collapse = " + ")))
formula.soli

# fit full hurdle models
hurdle.corrupt.out <- llply(data.mi.m9192, function(x)
  return(hurdle(formula.corrupt,
                dist = "negbin",
                zero.dist = "binomial", link = "logit",
                data = x)))
hurdle.soli.out <- llply(data.mi.m9192, function(x)
  return(hurdle(formula.soli,
                dist = "negbin",
                zero.dist = "binomial", link = "logit",
                data = x)))

# Get MI results
hurdle.corrupt.results <- combineHurdleMI(fitted.obj = hurdle.corrupt.out)
hurdle.corrupt.results
hurdle.soli.results <- combineHurdleMI(fitted.obj = hurdle.soli.out)
hurdle.soli.results

#######################################################
## Robustness Checks: Alternative Samples of Data
#######################################################
# fit full hurdle model (logit + negative binomial) to 10 MI datasets, period 1991-2010
hurdle.full.out <- llply(data.mi$imputations, function(x)
  return(hurdle(formula.full,
                dist = "negbin",
                zero.dist = "binomial", 
                link = "logit",
                data = x)))
hurdle.full.results <- combineHurdleMI(fitted.obj = hurdle.full.out)
hurdle.full.results

# negative binomial excluding adoption < 10 observations
nb.m10.full <- llply(data.mi.m10, function(x)
  return(glm.nb(formula.full,
                data = x,
                control = glm.control(epsilon = 1e-8, maxit = 200, trace = FALSE))))
nb.m10.full.results <- combineGlmMI(nb.m10.full)
nb.m10.full.results

#######################################################
## Fig 7. Negative Binomial Model of International Adoption, 
## Excluding Observations with Fewer than Ten Adoptions.
#######################################################
# set p-value
alpha <- 0.05

# extract coefficients
coef.gt10 <- nb.m10.full.results$coef

# set multiplier based on alpha and degrees of freedom
multiplier.gt10 <- qt(1 - alpha/2, df = nb.m10.full.results$dof)

# calculate high and low values
coef.gt10.high <- coef.gt10 + nb.m10.full.results$se*multiplier.gt10
coef.gt10.low <- coef.gt10 - nb.m10.full.results$se*multiplier.gt10

# omit factors
coef.gt10 <- coef.gt10[!grepl("factor", names(coef.gt10))]
coef.gt10.high <- coef.gt10.high[!grepl("factor", names(coef.gt10.high))]
coef.gt10.low <- coef.gt10.low[!grepl("factor", names(coef.gt10.low))]

# omit intercept
coef.gt10 <- coef.gt10[!grepl("Intercept", names(coef.gt10))]
coef.gt10.high <- coef.gt10.high[!grepl("Intercept", names(coef.gt10.high))]
coef.gt10.low <- coef.gt10.low[!grepl("Intercept", names(coef.gt10.low))]

# make list
coef.gt10.fig <- list("mean" = coef.gt10,
                      "high" = coef.gt10.high,
                      "low" = coef.gt10.low)

# assign variable names
coef.names.gt10 <- c("Requlatory Quality",
                     "Nationalism",
                     "ln(Cumulative Adoption)",
                     "Hague Adoption (both)",
                     "Language Commonality",
                     "ln(Migrant Stock)",
                     "Colonial Ties",
                     "ln(Distance)",
                     "ln(Youth Population)",
                     "ln(Real GDP per Capita)",
                     "Measles Immunization",
                     "Armed Conflict",
                     "Islamic Law")

# prepare parameters for plot
y.axis.coef.gt10 <- c(length(coef.gt10.fig$mean):1)

# open pdf device
pdf("coefPlotGT10.pdf", height = 5, width = 7) 
par(mar = c(2, 9, 2, 1))
plot(coef.gt10.fig$mean, 
     y.axis.coef.gt10, 
     type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = .6, col = "black", 
     xlim = c(round(min(coef.gt10.fig$low), digits = 2), 
              round(max(coef.gt10.fig$high), digits = 2) + 0.06),
     ylim = c(min(y.axis.coef.gt10), 
              max(y.axis.coef.gt10)),
     xaxs = "r", main = "") 
segments(coef.gt10.fig$low, 
         y.axis.coef.gt10, 
         coef.gt10.fig$high, 
         y.axis.coef.gt10, 
         lwd =  1.5, 
         col = "black")
axis(1, at = round(seq(round(min(coef.gt10.fig$low), digits = 2),
                       round(max(coef.gt10.fig$high), digits = 2) + 0.06,
                       by = 0.1), digits = 2),
     labels = round(seq(round(min(coef.gt10.fig$low), digits = 2),
                        round(max(coef.gt10.fig$high), digits = 2) + 0.06,
                        by = 0.1), digits = 2),
     tick = TRUE, 
     cex.axis = .8, 
     mgp = c(2, .5, 0)) 
axis(2, at = y.axis.coef.gt10, 
     label = coef.names.gt10, 
     las = 1, 
     tick = TRUE,
     cex.axis = .8) 
segments(0, 0, 0, 17, lty = 2)
mtext("Excluding < 10 Adoption Obs.", side = 3)
text(0.3, 3, paste("n = ", nb.m10.full.results$n), cex = 0.8)
text(0.379, 2, paste("logLik = ", round(nb.m10.full.results$loglik, digits = 2)), cex = 0.8)
text(0.35, 1, paste("AIC = ", round(nb.m10.full.results$aic, digits = 2)), cex = 0.8)
dev.off()

#######################################################
## Table B.2. Descriptive Statistics.
#######################################################
# select variables included in descriptive statistics
# use the subset of data, 1993-2010
des.vars <- list(data.9310$adopt_tot,
                 data.9310$reg_quality_o_L1,
                 data.9310$control_corrupt_o_L1,
                 data.9310$execnat_o_L1,
                 data.9310$Mjussolis_o_L1,
                 data.9310$adopt_o_cumsum_L1,
                 data.9310$hague_adopt_entry_both_L1,
                 data.9310$comlang_ethno,
                 data.9310$Istock_L1/1000,
                 data.9310$colony,
                 data.9310$distcap,
                 data.9310$pop14_mil_L1,
                 data.9310$rgdppc_exp_L1/1000,
                 data.9310$immune_measles_L1,
                 data.9310$n_mortality_L1,
                 data.9310$conflict_o_L1,
                 data.9310$islam_law_o)

# assign full variable names
des.nam <- c("Adoption", "Regulatory Quality", "Control of Corruption",
             "Nationalism", "Jus Soli",
             "Cumulative Adoption", "Hague Adoption (both)",
             "Language Commonality", "Migrant Stock (thousands)", "Colonial Ties", "Distance (km)",
             "Youth Population", "Real GDP per Capita (thousands)", "Measles Immunization", "Infant Mortality", "Armed Conflict", "Islamic Law")

# output descriptive stats table
des.tab <-  tableContinuous(vars = des.vars,
                            nams = des.nam,
                            stats = c("mean", "min", "max", "n", "na"),
                            prec = 2,
                            col.tit.font = c("bf", "", "sf", "it", "rm"),
                            cap = "Descriptive Statistics", lab = "tab:des",
                            font.size = "normalsize", longtable = FALSE)


#######################################################
## Fig B.1. Covariate Correlation Matrix.
#######################################################
library(polycor)
library(reshape2)
library(mgcv)

# create vectors of covariates
ivars.polycor <- c("reg_quality_o_L1", "control_corrupt_o_L1",
                   "execnat_o_L1", "Mjussolis_o_L1", "adopt_o_cumsum_L1", "hague_adopt_entry_both_L1",
                   "comlang_ethno", "Istock_L1", "colony", "distcap",
                   "pop14_mil_L1", "rgdppc_exp_L1", "immune_measles_L1",
                   "conflict_o_L1", "islam_law_o")

# subset
plot.df <- data.mi.m9192[[1]][, colnames(data.mi.m9192[[1]]) %in% ivars.polycor]

# rename
names(plot.df) <- c("Cumulative Adoption", "Language Commonality", "Colonial Ties", "Distance", "Migrant Stock",
                    "Armed Conflict",
                    "Regulatory Quality", "Control of Corruption", "Real GDP per Capita",
                    "Nationalism", "Jus soli", "Youth Population", "Measles Immunization", "Islamic Law", "Hague Adoption (both)")

# plot
# reshape to long format
plot.df <- melt(hetcor(plot.df, use = "pairwise.complete.obs")$correlations)

# reorder correlation order
plot.df <- plot.df[order(plot.df$value), ]

# reorder variable order
plot.df$Var1 <- reorder(plot.df$Var1, plot.df$value)
plot.df$Var2 <- reorder(plot.df$Var2, plot.df$value)

# plot heatmap
p <- ggplot(data = plot.df, aes(x = Var1, y = Var2, fill = value))
p <- p + geom_tile()
p <- p + scale_fill_gradient2(name = "Correlation", breaks = seq(-1, 1, by = .25),
                              space = "Lab")
p <- p + guides(fill = guide_colorbar(barwidth = .75, ticks = FALSE))
p <- p + labs(x = NULL, y = NULL)
p <- p + theme(axis.text.x = element_text(angle = 90, hjust = 1))
pdf(file = "corFig.pdf", width = 8, height = 8)
p
dev.off()

###########################################################################
## Fig B.2. Excess Zeros and Clustering in Directed Dyad Adoption Data.
###########################################################################
# plot histogram for dyadic adoption counts
# use the subset of data, 1993-2010
adopt.hist <- ggplot(data.9310, aes(x = adopt_tot)) +
  geom_histogram(binwidth = 1) +
  #scale_x_sqrt(breaks = c(0, 1, 1000, 2000, 4000, 6000, 8000)) +
  xlab("Total Adoption (omitting > 100 for presentation purpose)") + ylab("Count") +
  coord_cartesian(xlim = c(0, 10^2)) +
  theme_bw()

# plot boxplot of total adoptions by destination country
adopt.bp.d <- ggplot(data.9310, aes(factor(iso3c_d), asinh(adopt_tot))) +
  geom_boxplot() +
  xlab("Country") + ylab("Total Adoption \n (inverse hyperbolic sine transformation)") +
  theme_bw()

# save
pdf(file = "adopt.pattern.9310.pdf", width = 7, height = 7)
grid.arrange(adopt.hist, adopt.bp.d)
dev.off()


####################################################################################################
## Table D.1. Fitted Hurdle Negative Binomial Model Results with Multiple Im- putation, 1993-2010.
####################################################################################################
# output results in LaTeX table format
texreg(list(hurdle.full.m9192.results$texreg.zero, hurdle.full.m9192.results$texreg.count),
       custom.model.names = c("Hurdle:$Prob(Y>0)$","Hurdle:$E(Y)$"),
       custom.coef.names = c("Intercept",
                             "Regulatory Quality",
                             "Nationalism",
                             "ln(Cumulative Adoption)",
                             "Hague Adoption (both)",
                             "Language Commonality",
                             "ln(Migrant Stock)",
                             "Colonial Ties",
                             "ln(Distance)",
                             "ln(Youth Population)",
                             "ln(Real GDP per Capita)",
                             "Measles Immunization",
                             #"Infant Mortality",
                             "Armed Conflict",
                             "Islamic Law",
                             "factor(ifs_d)112", "factor(ifs_d)124", "factor(ifs_d)128", "factor(ifs_d)132", "factor(ifs_d)134", "factor(ifs_d)136",
                             "factor(ifs_d)138", "factor(ifs_d)142", "factor(ifs_d)144", "factor(ifs_d)146", "factor(ifs_d)156", "factor(ifs_d)172",
                             "factor(ifs_d)176", "factor(ifs_d)178", "factor(ifs_d)184", "factor(ifs_d)193", "factor(ifs_d)196", "factor(ifs_d)436",
                             "factor(year)1994", "factor(year)1995", "factor(year)1996", "factor(year)1997",
                             "factor(year)1998", "factor(year)1999", "factor(year)2000", "factor(year)2001", "factor(year)2002", "factor(year)2003",
                             "factor(year)2004", "factor(year)2005", "factor(year)2006", "factor(year)2007", "factor(year)2008", "factor(year)2009",
                             "factor(year)2010"),
       stars = c(0.001, 0.01, 0.05, 0.1), symbol = "\\wedge", digits = 3, omit.coef = "factor",
       #scriptsize=TRUE, 
       single.row = FALSE, float.pos = "!ht",
       caption = "Fitted Hurdle Negative Binomial Model Results with Multiple Imputation, 1993-2010.",
       label = "tb:fm_hurdle_full", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)

####################################################################################################
## Table D.2. Substantive Effect Estimates with Simulated 95% Confidence Intervals
####################################################################################################
# combine previously simulated effects in a dataframe
effects.table <- data.frame(effect.names,
                            effect.zero.vec, effect.zero.low.vec, effect.zero.high.vec,
                            effect.negbin.vec, effect.negbin.low.vec, effect.negbin.high.vec)

# assign names
names(effects.table) <- c("Covariates", "Prob. Estimates", "2.5%", "97.5%", "Count Estimates", "2.5%", "97.5%")

# convert to xtable format
effects.xtable <- xtable(effects.table, 
                         caption = "Substantive Effect Estimates with Simulated 95\\% Confidence Intervals", 
                         label = "tb:first_difference", digits = 3)

# print table in LaTeX format
print.xtable(effects.xtable, floating = TRUE, table.placement = "!ht", 
             caption.placement = "bottom", include.rownames = FALSE, include.colnames = TRUE, booktabs = TRUE)


####################################################################################################
## Table E.1. Robustness Models 1: Poisson, Negative Binomial, Zero-Inflated Negative Binomial Models.
####################################################################################################
# output results in LaTeX table format
texreg(list(pois.full.results$texreg,
            nb.full.results$texreg,
            zi.full.results$texreg.zero, zi.full.results$texreg.count),
       custom.model.names = c("Poisson", "NegBin", "ZINB:$Prob(Y=0)$","ZINB:$E(Y)$"),
       custom.coef.names = c("Intercept",
                             "Regulatory Quality",
                             "Nationalism",
                             "ln(Cumulative Adoption)",
                             "Hague Adoption (both)",
                             "Language Commonality",
                             "ln(Migrant Stock)",
                             "Colonial Ties",
                             "ln(Distance)",
                             "ln(Youth Population)",
                             "ln(Real GDP per Capita)",
                             "Measles Immunization",
                             #"Infant Mortality",
                             "Armed Conflict",
                             "Islamic Law",
                             "factor(ifs_d)112", "factor(ifs_d)124", "factor(ifs_d)128", "factor(ifs_d)132", "factor(ifs_d)134", "factor(ifs_d)136",
                             "factor(ifs_d)138", "factor(ifs_d)142", "factor(ifs_d)144", "factor(ifs_d)146", "factor(ifs_d)156", "factor(ifs_d)172",
                             "factor(ifs_d)176", "factor(ifs_d)178", "factor(ifs_d)184", "factor(ifs_d)193", "factor(ifs_d)196", "factor(ifs_d)436",
                             "factor(year)1994", "factor(year)1995", "factor(year)1996", "factor(year)1997",
                             "factor(year)1998", "factor(year)1999", "factor(year)2000", "factor(year)2001", "factor(year)2002", "factor(year)2003",
                             "factor(year)2004", "factor(year)2005", "factor(year)2006", "factor(year)2007", "factor(year)2008", "factor(year)2009",
                             "factor(year)2010"),
       stars = c(0.001, 0.01, 0.05, 0.1), symbol = "\\wedge", digits = 3, omit.coef = "factor",
       #scriptsize=TRUE, 
       single.row = FALSE, float.pos = "!ht",
       caption = "Robustness Models 1: Poisson, Negative Binomial, Zero-Inflated Negative Binomial Models with Multiple Imputation",
       label = "tb:fm_robust1", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)

####################################################################################################
## Table E.2. Robustness Models 2: Fitted Hurdle Negative Binomial Models with Alternative Indicators.
####################################################################################################
# output results in LaTeX table format
texreg(list(hurdle.corrupt.results$texreg.zero, hurdle.corrupt.results$texreg.count,
            #hurdle.infant.results$texreg.zero, hurdle.infant.results$texreg.count,
            hurdle.soli.results$texreg.zero, hurdle.soli.results$texreg.count),
       custom.model.names = c("Hurdle:$Prob(Y>0)$","Hurdle:$E(Y)$",
                              #"Hurdle:$Prob(Y>0)$","Hurdle:$E(Y)$",
                              "Hurdle:$Prob(Y>0)$","Hurdle:$E(Y)$"),
       custom.coef.names = c("Intercept",
                             "Control of Curruption",
                             "Nationalism",
                             "ln(Cumulative Adoption)",
                             "Hague Adoption (both)",
                             "Language Commonality",
                             "ln(Migrant Stock)",
                             "Colonial Ties",
                             "ln(Distance)",
                             "ln(Youth Population)",
                             "ln(Real GDP per Capita)",
                             "Measles Immunization",
                             "Armed Conflict",
                             "Islamic Law",
                             "factor(ifs_d)112", "factor(ifs_d)124", "factor(ifs_d)128", "factor(ifs_d)132", "factor(ifs_d)134", "factor(ifs_d)136",
                             "factor(ifs_d)138", "factor(ifs_d)142", "factor(ifs_d)144", "factor(ifs_d)146", "factor(ifs_d)156", "factor(ifs_d)172",
                             "factor(ifs_d)176", "factor(ifs_d)178", "factor(ifs_d)184", "factor(ifs_d)193", "factor(ifs_d)196", "factor(ifs_d)436",
                             "factor(year)1994", "factor(year)1995", "factor(year)1996", "factor(year)1997",
                             "factor(year)1998", "factor(year)1999", "factor(year)2000", "factor(year)2001", "factor(year)2002", "factor(year)2003",
                             "factor(year)2004", "factor(year)2005", "factor(year)2006", "factor(year)2007", "factor(year)2008", "factor(year)2009",
                             "factor(year)2010",
                             "Regulatory Quality",
                             #"Infant Mortality",
                             "Jus soli"),
       stars = c(0.001, 0.01, 0.05, 0.1), symbol = "\\wedge", digits = 3, omit.coef = "factor",
       #scriptsize=TRUE, 
       single.row = FALSE, float.pos = "!ht",
       caption = "Robustness Models 2: Fitted Hurdle Negative Binomial Models with Alternative Measures",
       label = "tb:fm_robust2", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)

####################################################################################################
## Table E.3. Robustness Models 3: Fitted Hurdle Negative Binomial Models with Alternative Samples of Data
####################################################################################################
# output results in LaTeX table format
texreg(list(hurdle.full.results$texreg.zero, hurdle.full.results$texreg.count,
            nb.m10.full.results$texreg),
       custom.model.names = c("Hurdle:$Prob(Y>0)$","Hurdle:$E(Y)$",
                              "NegBin"),
       custom.coef.names = c("Intercept",
                             "Regulatory Quality",
                             "Nationalism",
                             "ln(Cumulative Adoption)",
                             "Hague Adoption (both)",
                             "Language Commonality",
                             "ln(Migrant Stock)",
                             "Colonial Ties",
                             "ln(Distance)",
                             "ln(Youth Population)",
                             "ln(Real GDP per Capita)",
                             "Measles Immunization",
                             #"Infant Mortality",
                             "Armed Conflict",
                             "Islamic Law",
                             "factor(ifs_d)112", "factor(ifs_d)124", "factor(ifs_d)128", "factor(ifs_d)132", "factor(ifs_d)134", "factor(ifs_d)136",
                             "factor(ifs_d)138", "factor(ifs_d)142", "factor(ifs_d)144", "factor(ifs_d)146", "factor(ifs_d)156", "factor(ifs_d)172",
                             "factor(ifs_d)176", "factor(ifs_d)178", "factor(ifs_d)184", "factor(ifs_d)193", "factor(ifs_d)196", "factor(ifs_d)436",
                             "factor(year)1992", "factor(year)1993", "factor(year)1994", "factor(year)1995", "factor(year)1996", "factor(year)1997",
                             "factor(year)1998", "factor(year)1999", "factor(year)2000", "factor(year)2001", "factor(year)2002", "factor(year)2003",
                             "factor(year)2004", "factor(year)2005", "factor(year)2006", "factor(year)2007", "factor(year)2008", "factor(year)2009",
                             "factor(year)2010"),
       stars = c(0.001, 0.01, 0.05, 0.1), symbol = "\\wedge", digits = 3, omit.coef = "factor",
       #scriptsize=TRUE, 
       single.row = FALSE, float.pos = "!ht",
       caption = "Robustness Models 3: Fitted Hurdle Negative Binomial Models with Alternative Samples of Data",
       label = "tb:fm_robust3", use.packages = FALSE, booktabs = TRUE, dcolumn = TRUE)
